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ABSTRACT 

We present a new three-dimensional radiative transfer (RT) code, RADAMESH, 
based on a ray-tracing, photon-conserving and adaptive (in space and time) scheme. 
RADAMESH uses a novel Monte Carlo approach to sample the radiation field within the 
computational domain on a "cell-by-cell" basis. Thanks to this algorithm, the compu- 
tational efforts are now focused where actually needed, i.e. within the Ionization-fronts 
(I- fronts). This results in an increased accuracy level and, at the same time, a huge gain 
in computational speed with respect to a "classical" Monte Carlo RT, especially when 
combined with an Adaptive Mesh Refinement (AMR) scheme. Among several new 
features, RADAMESH is able to adaptively refine the computational mesh in correspon- 
dence of the I-fronts, allowing to fully resolve them within large, cosmological boxes. 
We follow the propagation of ionizing radiation from an arbitrary number of sources 
and from the recombination radiation produced by H and He. The chemical state of 
six species (HI, HII, Hel, Hell, Helll, e) and gas temperatures are computed with a 
time-dependent, non-equilibrium chemistry solver. We present several validating tests 
of the code, including the standard tests from the RT Code Comparison Project and 
a new set of tests aimed at substantiating the new characteristics of RADAMESH. Using 
our AMR scheme, we show that properly resolving the I-front of a bright quasar dur- 
ing Reionization produces a large increase of the predicted gas temperature within the 
whole HII region. Also, we discuss how H and He recombination radiation is able to 
substantially change the ionization state of both species (for the classical Stromgren 
sphere test) with respect to the widely used "on-the-spot" approximation. 

Key words: 

radiative transfer - methods: numerical - HII regions - intergalactic medium - diffuse 
radiation - cosmology:theory 



1 INTRODUCTION 

The interaction between matter and radiation is one of the 
fundamental mechanisms shaping the distribution of the 
baryonic component in the Universe, from stellar to cos- 
mological scales. This process often couples wildly different 
scales and the large dynamical range makes accurate mod- 
elling difficult. A notable example is given by the reioniza- 
tion of cosmic hydrogen at redshift 6 < z < 12 (for a re- 
view, see, e.g., Meiksin 2009). During this epoch, the HII 
regions generated by the first ionizing sources expand to in- 
tergalactic scales and overlap, leaving most of the Universe 
highly ionized and re-heated by several thousand Kelvin de- 
grees. The brightest, high-redshift quasars, may have pro- 
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duced even larger HII regions before the end of Reionization, 
with linear sizes extending up to a hundred comoving Mpc. 
Nonetheless, most of the evolution in the gas temperature 
and ionization state still happened on much smaller scales, 
i.e. within the Ionization-fronts, whose sizes are comparable 
to the local photon mean-free-path just outside the HII re- 
gions, about three orders of magnitude smaller than the HII 
regions themselves for gas at mean cosmic density. 

Even without considering such an extreme case in dy- 
namical range, the numerical solution of the full radiative 
transfer equations still represents a big computational chal- 
lenge. This is mostly due to two factors. The first is the 
high-dimensionality: seven variables are required for a full 
specification of the radiation field (three spatial variables, 
two angular directions, photon frequency, and time). The 
second is the intrinsic "non-locality" of the problem: the 
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radiation field at a given point is determined by the gas 
properties along the lines of sight towards all the sources of 
radiation (including the diffuse medium itself). 

For these reasons, current numerical models rely on ap- 
proximations aimed to decrease the problem dimensionality 
and, at the same time, the computational costs associated 
with the "non-locality" of the full radiation transfer (e.g., 
Abel, Norman & Madau 1999, Gnedin & Abel 2001, Maselli, 
Ferrara & Ciardi 2003, Razoumov & Cardall 2005, Mellema 
et al. 2006, Rijkhorst et al. 2006, Ritzerveld & Icke 2006, 
Whalen & Norman 2006; for a review, see also Iliev et al. 
2006 and references therein; more recent codes include, e.g., 
Trac & Cen 2007, Semelin et al. 2007, Aubert & Teyssier 
2008, Altay et al. 2008, Pawlik & Schaye 2008, Finlator et 
al. 2009, Petkova & Springel 2009). A widely used approx- 
imation is the so called "method of characteristics" which 
is particularly suited for cases where the radiation field is 
completely dominated by individual sources (with no con- 
tribution from the radiative recombinations produced by 
the medium itself). In this scheme, the radiation field is 
determined simply from the column densities along a line 
of sight between the sources and a given point in the com- 
putational domain, discretized in cells (or particles). In the 
"long-characteristics" flavour, the column densities are di- 
rectly calculated with a single ray from the sources to the 
cell centre, summing up all the contributions cut through 
the cells in between. In the faster, but less accurate, "short- 
characteristics" version, the column densities are determined 
in an ordered way starting from the cells closer to the source 
and moving outwards after summing up (using some inter- 
polation scheme) previous contributions. In both cases, sam- 
pling the radiation field with a single ray from the sources to 
the cell centre may be problematic in particular situations, 
i.e. when the cells are optically thick and in presence of 
strong density gradients. Practically, the consequences may 
be the loss of photon conservation and numerical artifacts 
in the I-front shapes (see, e.g., Mellema et al. 2006). 

In a cosmological situation, e.g. for the reionization sce- 
nario described before, we would like to have an accurate 
method for the very typical case in which numerical reso- 
lution imposes optically thick cells. This may be obtained 
sampling the radiation field in a statistical, homogeneous 
way, using the widely used Monte Carlo techniques. In this 
case, a series of rays is randomly casted around the sources 
in order to obtain the right solid angle distribution, averag- 
ing the contributions to the radiation field within individ- 
ual (three-dimensional) volume elements. The drawbacks of 
this method are the computational costs: many rays (per 
cell) are required to avoid statistical noise and the method 
becomes less and less efficient as one moves outwards from 
the sources (but see, Abel & Wandelt 2002). Moreover, if 
the grid is composed by cells of different sizes, the required 
number of rays to be casted is determined by the size of 
the smallest cells, irrespectively of the volume fraction they 
actually occupy. 

In this work, we develop a novel approach aimed to ob- 
tain an accurate solution of the radiative transfer problem 
in cosmological situations in an efficient way: combining the 
Monte Carlo scheme with the "cell-by-cell" approach typi- 
cal of the characteristic methods. As we mentioned at the 
beginning of this section, despite the large scales associated 
with, e.g., the reionization process, the medium properties 



are very often determined within a small fraction of its vol- 
ume, i.e. within the I-fronts. The basic idea is to concen- 
trate the computational efforts onto this part of the com- 
putational volume, where they are actually needed. This is 
achieved thanks to a new algorithm which casts, with the 
correct solid angle distribution, a series of rays for each cell 
individually. In this way, we are able to use an algorithm 
that adaptively determines: i) the number of rays needed 
for a particular cell to achieve the convergence of the radi- 
ation field (irrespectively of the cell size, or distance from 
the source), ii) in which part of the volume to apply the 
(expensive) Monte Carlo ray-tracing. The result is that the 
computing time of the RT scales now with the number of 
cells to be evolved during a simulation-step, i.e., typically 
the cells contained within the I-fronts, with a huge gain in 
computational speed (and accuracy) with respect to a clas- 
sical Monte Carlo. 

With the new "cell-by-cell" approach we are also able 
to gain the full advantages of Adaptive Mesh Refinement 
(AMR) methods (see, e.g. Berger & Oliger 1984) that in- 
crease the spatial resolution where needed. In particular, we 
implemented a scheme, based on the original clustering al- 
gorithm by Berger & Rigoutsos (1991), to adaptively refine 
the mesh in correspondence of the I-fronts during the course 
of the simulation. This allows us to fully resolve the small 
scales associated with the I-fronts within the large, cosmo- 
logical simulation boxes that follow the reionization process. 
Resolving the I-fronts has a profound impacts for a large 
range of astrophysical applications, e.g., to accurately pre- 
dict the temperature state of the gas surrounding a bright 
quasar during Reionization. In the validating test section 
(section 4), we will present an example of how resolving the 
I-front may be important in the determination of the gas 
temperature within the whole HII region. 

The paper is organized as follows. In section 2, we 
review the basic radiative transfer equation. The com- 
putational method used by our new RT code, RADAMESH 
(Radiative-transfer on ADAptive MESH), is presented in 
section 3, In section 4, we show the validating tests of the 
code. We conclude in section 5. 



2 BASIC EQUATIONS 

The cosmological radiative transfer equation in comoving 
coordinates (e.g., Norman, Paschos & Abel 1998) is given 
by 

idi„ , h-vi u H(t) ar„ ... 

-^t-H z [y-= 3i„) = k„(S v -lu) , (1) 

c at a c ov 

where I v = I(t, x, h, v) is the monochromatic specific in- 
tensity of the radiation field, n is a unit vector along the 
direction of propagation of the ray, H(t) = a/a is the (time- 
dependent) Hubble constant, a = 1 ^° m is the ratio of cos- 
mic scale factors between photon emission at frequency v 
and the present time t, k v denotes the opacity at frequency 
v and S v is the source function of the medium. 

If the scale of interest L is much smaller than the Hub- 
ble radius, c/H, (as always in our case) and the medium 
properties are changing on a time scale shorter than the 
light crossing time L/c, equation |TJ reduces to the classi- 
cal, static radiative transfer equation: 
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n-VI v = k v {S v - I v ) . 
This equation admits the following solution, 



(2) 



I v {j v ,n) = I„(p,n) ■ e r " + / S v (r v ,n)dT' v , (3) 

Jo 

where r v is the optical depth along n. This solution is in gen- 
eral adequate for cosmological simulations except on small 
distances from the sources (or, equivalently, for very bright 
sources). In this case, the approximations break down allow- 
ing the Ionization-fronts to expand faster than the speed of 
light. In section [3.3.1 1 we discuss an approximate correction 
to solve this unphysical behaviour. It is often convenient to 
express I v as a sum of the attenuated, direct radiation from 
individual sources (I„ 1T ) and the diffuse radiation (J^ 1 ) gen- 
erated within the medium: 



l v {T v ,h) = lt"(r v ,h) +/„' (r„,n) 



(4) 



Using equation the radiation field intensity can be 
specified at any point of the simulation box given the value 
of the optical depth to the individual sources and the source 
term of the medium. From the radiation field, it is then 
possible to derive for each species i the photoionization rate 
(per particle): 



r 4 (5) = 



1 ») dl/dn 

hv 



(5) 



and the gas photoheating rate (per unit volume) 
d (a?) = n 



aiiv) 1 "^'™' h{v - Vi)dvdQ , (6) 

where Vi and Oi are the frequency threshold and the 
cross-section for the ionization of species i, respectively, and 
i%i is the physical number density. Analogously to eq. Q, 
we can split Fi and Gi into the sum of direct (Tf 1T ,Gf lr ) and 
diffuse (rf ft ,Gf ft ) components. 

Finally, Ti and Gi are used to compute the chemistry 
evolution of the neutral fraction of hydrogen (/hi), neutral 
helium (/hci), singly ionized helium (/hcIi), and the tem- 
perature (expressed in terms of the total energy density 
E = (3/2)ntotkbT) according to the following rate equa- 
tions: 



dfm = 
dt 

d/Hd 

dt 

d/HoII 

dt 

dE 
lit 



-/hiThi — n c /Hi/3fflC + n c fmiamC 



(7) 



-/HelTHel — We/Hcl/^HcI C + 7l c /Hell (tVHcI + £)C 

(8) 

— /HcIirHcII — "c/h c Ii/3hcIiC + ne/HcIIILtHcIlC 



+ 2HS = ^]n l /,G I -A(n I ,T) 



(9) 
(10) 



Here cti, Pi and £ are the temperature-dependent radia- 
tive recombination (to all levels), collisional ionization, and 
dielectronic recombination coefficients, respectively; n c — 
"-h(1 — /hi) + "-ho (2 — 2/hcI — /hcIi) is the number density of 
electrons, C = (n 2 ) / (n) 2 is the gas clumping factoiQ , and 

1 For simplicity, we assume the same value of C for each species 
in the current version of the code. 



A is the total cooling rate (including recombinations, col- 
lisional excitations, Compton, brehmstrahlung and Hubble 
cooling). We use the analytical fits of Hui & Gnedin (1998) 
for the ionization, recombination and cooling rates. 

One of the main challenges in the calculation of the pho- 
toionization/heating rate is represented by the presence of 
the diffuse term in eq. Q, since it requires the transport of 
the diffuse photons generated from atomic recombinations 
in every point of the medium. A widely used solution to 
this problem is the so called "on-the-spot" approximation 
(OTS or Case B, originally proposed by Baker & Menzel 
1938), which assumes that every diffuse, ionizing photon is 
absorbed exactly in the same point where it is generated 
(i.e., the photon mean-free-path is much smaller than the 
resolution scale). In this case, assuming for the moment a 
pure hydrogen medium, the diffuse part of the photoioniza- 
tion rate can be written as: 



Fhi = n e /HiaHiC = n e /Hi(am — Qhi)C 



(11) 

where a HI is the recombination coefficient to the hydrogen 
ground level, and a H i is the analogous coefficient for the 
levels I > 1. Substituting this relation in eq. Q, we obtain: 

-jp = — /kiThY — n c fnif3mC + n e fmiamC , (12) 

i.e., a relation that depends only on the photoionization 
rate from individual sources (FhY) and, thus, greatly sim- 
plifies the radiative transfer problem. Analogous relations 
can be obtained for Hel and Hell, assuming that every ion- 
izing photon from recombinations is absorbed by the same 
species from which it originates (obviously, this approxima- 
tion breaks down if we consider that non-self-ionizing pho- 
tons from Hell recombinations can actually ionize both Hel 
and HI, see section 4.8 for a detailed discussion). 

As also noticed by Ritzerveld (2005), despite being 
widely used, the OTS approximation is based on an incor- 
rect argument: if the local mean- free-path is very small as 
assumed, i.e., the local optical depth is very high for the dif- 
fuse photons generated locally, even higher the optical depth 
would be for the radiation coming from discrete, non-local 
sources (unless their spectral energy distribution is much 
harder). For instance, in the classical Stromgren sphere sit- 
uation with a monochromatic source (see Test 1 below), the 
regime where the OTS approximation may hold corresponds 
to regions where the directional flux from the central source 
is not able to penetrate, i.e. only at the extreme edge of the 
Stromgren sphere. The opposite case (called Case A) is rep- 
resented by the situation in which the mean-free-path of the 
diffuse photons goes to infinity and, effectively, F^i — 0. 
Despite the loss of photon conservation for a bounded re- 
gion like a Stromgren sphere, this approximation is actually 
more correct than the OTS in most situations as we will 
show later. 



3 COMPUTATIONAL METHOD 

The computational volume in RADAMESH is discretized in a 
(series of) regular, block-structured grid(s) where the phys- 
ical properties of the medium are associated with a zone- 
centered grid element, i.e. with a cell. The current possible 
choices for the computational domain in RADAMESH are: i) 
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single, uniform Cartesian grid, ii) static multi-mesh struc- 
ture, iii) evolving multi-mesh structure. In the static case, 
the grid is spatially fixed to the initial, single or multi-mesh 
structure obtained, e.g., from the output of an Adaptive 
Mesh Refinement (AMR) hydro-simulation. In the evolving 
case, this initial structure is adaptively refined (unrefined) 
each time the cells satisfy a chosen refinement (unrefine- 
ment) criterion. The refining procedure is described in de- 



tails in section 3.3.4 The radiation field (from individual 
sources and from diffuse radiation) is discretized in a series of 
rays that are propagated through this single or multi-mesh 
structure using a ray-tracing algorithm. Before discussing 
the ray-tracing method, it is necessary to describe more in 
detail the multi-mesh implementation in RADAMESH. 

3.1 Multi-mesh in RADAMESH 

Multi-mesh domains in RADAMESH are composed of a nested 
hierarchy of rectangular grids of different sizes and levels 
of refinement, following the implementation called "patch- 
based AMR", originally described in Berger & Oliger (1984). 
The position and aspect ratios of this patches are optimized 
in order to increase the spatial (and temporal) resolution 
where needed within the computational box. Other possi- 
ble multi-mesh implementations have been proposed in the 
past, like, e.g, the "tree-based AMR" (see Khokhlov 1998 
and references therein). In this case, each cell (or n 3 group 
of cells in three-dimensions; typically n = 2) is refined into 
children cells, on a cell-by-cell basis, and a "grid" is the 
analogous of a single cell (or a n 3 group of cells). These 
methods were originally developed for the numerical solu- 
tion of hydrodynamical equation and they have advantages 
and disadvantages with respect to each other for this par- 
ticular problem. In our case, patch-based AMR presents a 
significant advantage with respect to other implementations: 
it reduces the total number of grids. One of the major prob- 
lem presented by a ray-tracing algorithm is the fact that a 
significant part of the computational time may be spent in 
"crossing" the simulation box. This is in a sense unavoid- 
able, because solving the radiative transfer equations is a 
highly non-local problem. Optimizing the data and memory 
structure is thus fundamental in order to have an efficient 
algorithm. Minimizing the number of "crossing" grids helps 
reducing the overhead associated with the ray-tracing in a 
multi-mesh structure. 

In order to increase the efficiency in "crossing" the com- 
putational volume, grids and cells in RADAMESH at different 
levels of refinement are connected with two separate hierar- 
chical trees and linked lists. Grids are strictly nested, i.e., 
each grid at level I is fully contained by a single grid at level 
/ — 1 (the "parent" grid). Note that, although this may in- 
crease the total number of grids, it allows a more efficient 
tree-search for both grids and cells. Grids that share the 
same "parent" are connected via a linked list. Thus, each 
grid gi at level I may have a maximum of three associated 
pointers: i) the "parent" grid at level I — 1, ii) the "next 
brother" grid in the linked list at level I, iii) the "son" grid 
at level I + 1 (i.e., the head of the linked list of grids fully 
contained within gi). Grids without an associated "son" are 
called "leaf" grids. This hierarchical tree is a light structure 
that requires small amounts of memory (if the number of 
grids is significantly smaller than the total number of cells, 



as it is always the case for "patch-based AMR") and that 
can be reconstructed or saved easily given its strict nesting. 
If refined, a cell at level I contains a group of rf subcells, 
where r; is the refinement factor. In principle, n may vary 
with levels, although commonly is fixed to a constant value 
(typically rj = 2). Analogously to the grid case, a cell with- 
out refinement is called a "leaf" cell. 

Combined with the grid-tree, there is another simple 
hierarchical structure in RADAMESH associated directly with 
the cells: each refined cell at level I points to the grid at 
I + 1 that contains its rf subcells. This simple (but more 
memory-consuming) cell-tree is built efficiently when needed 
from the grid-tree and, for this reason, is not saved during 
output or back-up in order to save memory usage. The use 
of the cell-tree allows to recover very quickly the location of 
the relevant subgrid (and thus the subcell) at level l + l when 
crossing the tree towards higher levels of refinement, without 
searching through the "next brother" linked list of the sub- 
grids. When searching for cells at coarser (or at the same) 
level, the grid-tree is used instead. The combination of the 
grid-tree and the simple cell-tree represents a good balance 
between flexibility (that translate into computational speed) 
of the box crossing algorithm and memory consumption of 
the tree-structure. 



3.2 Ray-tracing algorithm 

Ray-tracing is a well-studied problem in computational ge- 
ometry, with several efficient methods commonly used in 
computer graphics. One of the most efficient solution is the 
simple and fast grid-traversal algorithm of Amanatides & 
Woo (1987), particular suited in the case of a single, uniform 
Cartesian grid. Basically, this simple algorithm determines 
the intersection points between the boundaries of the cells 
and a ray (defined from its initial position and two direction 
angles) as it traverses the computational grid. 

We employ the grid and cell hierarchical trees discussed 
above to extend this fast algorithm from the single, uni- 
form grid to the multi-mesh structure. This is done simply 
adding to the original algorithm, after each cell boundary 
crossing, a (recursive) check for the need of changing the 
current traversing grid (either at the same or at a different 
level). Let us consider, for example, the case in which a ray 
is currently traversing a grid at level I = U. There are two 
situations in which a change of grid is required: a) the ray is 
still within the boundaries of the current grid but is entering 
a subgrid at level I > U ; b) the ray is exiting the boundaries 
of the current grid. We know that we are in the first situation 
simply checking the current cell tree: if the cell is associated 
with a subgrid, we immediately recover from the tree the 
memory location of the new grid at I = h + 1 and, from 
the intersection coordinates, the corresponding cell within 
this grid. We recursively continue the search until we reach 
a "leaf" cell. In the second case, the search consists of two 
phases. First, we cross the grid tree towards the "parent" 
grids at level I < U until we find a grid where the ray is 
fully contained (and not at the grid boundaries) or where 
the ray is entering the grid boundaries. Note that usually 



2 this is simply determined from the ray direction. 
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this search is very fast given that, for a typical nested struc- 
ture, we only need to go to the level I = U — 1. When the 
grid at I < U is found, we check if the cell where the ray is 
currently located is a "leaf". If this is not the case, we apply 
the same procedure as in case a). 

When traversing a "leaf" grid we can apply the original 
traversal algorithm without the need for the recursive check 
(and change of grids) described above. In this respect, the 
use of a "patch-based AMR", given its larger grids, result 
in a reduction of the computational costs of the traversing 
algorithm. 



3.3 RADAMESH step by step 

RADAMESH is based on a ray-tracing algorithm that propa- 
gates the (discretized) radiation field through the compu- 
tational volume. However, in the same spirit of the AMR 
method, the key element in RADAMESH is not represented by 
the ray or photon package but, instead, by the computa- 
tional cell. In other words, the propagation of the radiation 
field (and the solution of the radiative transfer equation) is 
performed on a "cell-by-cell" basis rather than "ray-by-ray" . 
Computationally, this translates into the change of the main 
loop in the algorithm, from rays to cells: the radiation field 
from the (discrete and diffuse) sources is propagated with 
a photon-conserving method separately for each cell in the 
computational volume. As we explain below, at the heart of 
the algorithm there is a new method that allows to draw a 
series of rays from a source, located at any position inside 
(or outside) the box, through a cubic cell with the correct 
solid angle distribution. This extends the "Monte Carlo" 
methods, where a series of rays is uniformly casted around 
a source to obtain the correct solid angle distribution, from 
the computational box to the cell level. 

This "cell-by-cell Monte Carlo" is a more natural ap- 
proach when most of the physical evolution of the medium 
happens rapidly only in a small portion of the simulated box 
at a time. This is a typical case in cosmological simulation, 
either because of a difference in density (e.g., self-shielded 
dense clumps in a mostly ionized medium) or in the ionized 
state (e.g., the expansion of an Ionization-front in a mostly 
neutral IGM during reionization) With the cell-based ap- 
proach, we can split the computational volume, either sin- 
gle or multi-mesh, grouping the cells with different evolving 
times or properties and focusing the computational efforts 
only where actually needed. This is obtained associating to 
each cell an individual time-step (At c ) that is used to deter- 
mine whether the cell can be considered active or not during 
the current simulation step. We discuss how we derive At c 
and how we define the active cell in section 13.3.21 

In detail, the computational algorithm consists of an 
iterative method divided into four main parts for each simu- 
lation step (and for each active cell): i) finding the ionization 
and photo-heating rates, ii) choosing the time-step, iii) solv- 
ing the chemistry equation for the evolution of the medium 
properties, and, if needed, iv) performing an adaptive re- 
finement of the computational mesh. The first part of the 
algorithm is the most time-consuming (and important) part 
of the method and is described in detail in the following 
section. 



3.3.1 Finding the ionization and photo-heating rates 

The photo-ionization (Ti) and photo-heating rates (Gi) are 
computed for each active, leaf cell (during the current time- 
step) with an iterative Monte Carlo procedure until conver- 
gence for both quantities is reached. The convergence level 
is typically set at 1%, but can be increased (decreased) de- 
pending on the requested level of accuracy. 

We compute separately, with two different methods, the 
contribution to the total photoionization/heating rate from 
discrete sources (rf lr ) and from diffuse emission (rf ff ).Most 
of the effort is dedicated to the calculation of (IV r ), since 
in many (cosmological) cases the directional flux from dis- 
crete sources represents the major contribution to the total 
photoionization rate. 

The procedure to obtain the values of Tf lr and Gf 1T is 
the following: for each iteration step, a packet of rays is prop- 
agated from selected points within the cell to the sources. 
These points are chosen with a Monte Carlo procedure based 
on a rejection algorithm that insures an uniform solid-angle 
distribution within the cell. The total flux is then re-scaled 
according to the fraction of the solid angle covered by the 
cell. We present in the Appendix the analytical approxima- 
tion that gives the solid angle of a (cubic) cell as seen by 
any point inside (or outside) the box. 

For each ray, we compute the HI, Hel, and Hell column 
densities both for the path length within the cell (AJVi) and 
for the path length connecting the cell edge to the source 
(Ni). The column densities are then used to calculate the 
frequency-dependent optical depths Ati(v) = ai(v)ANi and 
Ti{y) = o~i{v)Ni. Given the optical depths, we derive the 
probability that a photon with frequency v is absorbed by 
the species i within the cell: 

l-exp[- £ Ar»] 

Pity) = ~ 3 — {1 ~ expf-Ar^)]} . (13) 

£{l-exp[-Ar»]} 
j=i 

The optical depths and the spectral energy distribution of 
the sources are sampled into N v (logarithmically-spaced) 
frequency bins. The probability distributions Pi are used to 
compute the photoionization rate, including for the moment 
only the first term in eq. (J3j) , for each ray of the packet and 
for each source: 

ff 1 = J2 Pi^JiV^JexpI-nKJ] , (14) 

/,.=! 

where N 1 (y ilJ ) denotes the number of ionizing photons per 
unit time emitted by the source in the frequency bin i v . 
Finally, the value of r^ lr is obtained by averaging Ff' 1 over 
the number of rays in the packet, the cell volume and the 
fraction of the solid angle covered by the cell (see Appendix). 
Several packets of rays are generated until the value of r^ lr 
converges to the required level. The procedure is repeated 
for each source and the single values of r^ lr are added. The 
photoionization rates Gf 1T are obtained in a similar way, 
starting from 

£,dir = ^2 P x {v iv )h{v iv - i/ i )iV 7 (i/ il/ )expt-Tj(i'i„)] . (15) 

i„=l 

In order to avoid redundant calculations of the column 
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densities Ni from the cell edges to the sources for each ray 
packet, we first evaluate Ni on the cell vertexes "visible" 
from each source. If the difference between these values is 
less than a given (small) threshold, we do not use the full ray- 
casting algorithm. Instead, we use the Ni values at the ver- 
texes to interpolate (linearly) the needed values on the cell 
faces. Typically, the cells that need the full ray-casting algo- 
rithm are a few percent of the total volume and correspond 
to the Ionization-front regions (where the column-density 
varies rapidly). For the remaining cells, interpolating the 
column-densities is a good approximation and allows to sub- 
stantially speed-up the computation. Note that in this case, 
the algorithm is similar to a classical long-characteristic ray- 
tracing. 

To achieve a better performance it is also possible to 
choose a column-density threshold (or, equivalently, an op- 
tical depth threshold at the ionization limit) above which 
RADAMESH skips the calculation of the F dlr and G dlr values. 
This is computationally convenient, e.g., in the early stages 
of the expansion of an Ionization-front in a neutral and dense 
medium, when most of the volume is still optically thick to 
the source radiation. Analogously, it is possible to choose a 
minimum value of ANi and Ni below which an optically- 
thin approximation is used instead of the full ray-tracing 
algorithm. This is typically the situation for the ionized re- 
gion in proximity of the sources, where the optical depth is 
negligible. 

As mentioned in section 2, eq. Q permits I-fronts to 
propagate faster than light in near proximity of a bright 
source. An exact solution for this problem that does not 
require to fully solve the time-dependent radiative trans- 
fer equation and avoids the loss of photon conservation is 
only available in the case of an homogeneous medium (e.g., 
Shapiro et al. 2006) as we have also discussed elsewhere 
(Cantalupo et al. 2008). In the more general case of an in- 
homogeneous medium, we fix this unphysical behaviour by 
disregarding the ionizing flux of a given source at a distance 
r > (t s -c), where t s is the source current lifetime. Note that, 
the loss of photon conservation due to this approximation is 
typically restricted to the very early phases of the expansion 
of the I-front produced by very bright sources. 

Once the value of the directional photoioniza- 
tion/heating rate have been obtained, we calculate (with 
a simpler procedure) the diffuse component rates for the 
same, active cells. Also in this case we follow a "cell-by-cell" 
approach: given a particular (active) cell for which we want 
to compute r dlff and G dlff , we generate with a Monte Carlo 
procedure a set of rays that propagates outwards from the 
cell centre. For each cell encountered by the ray, we com- 
pute the medium opacity and the source term as a function 
of the n e , nun, riHell, n-Hein, and temperature. In particu- 
lar, we include: i) HI, Hel and Hell free-bound continuum 
(from Osterbrock 1989), ii) Hell Balmer continuum (from 
Ercolano & Storey 2006), iii) Hell two-photon continuum 
(from Nussbaumer & Schmutz 1984), and iv) the Hell Lya 
line. We neglect the Hel Balmer continuum and Hel emis- 
sion lines given their (relatively) small contribution to the 
total emissivity. Given the source function, we compute and 
add the respective r dlff and G dlff to the direct photoion- 
ization/heating rate, according to the equation presented in 
section 2. 

In this method, each ray represents a fixed fraction of 



the solid angle (0d) of the sky as seen by the cell, deter- 
mined by the chosen number of rays (that we leave in the 
current version as a free parameter). Note that this pro- 
cedure is accurate if the three-dimensional diffuse field is 
homogeneous over a scale corresponding to x r, where 
r is the distance from the cell. Obviously, if the medium 
is clumpy, this method becomes more and more inaccurate 
at increasing distances from the cell, unless the number of 
rays is increased accordingly. Note, however, that in most 
cosmological situations, the diffuse field represents a small 
component with respect to the direct radiation from sources. 

3.3.2 Choosing the time-step and the active cells 

As mentioned above, we associate an individual time-step 
At c to each cell in the computational volume. We first esti- 
mate At c for the cells with the newly calculated Fi by taking 
a fixed fraction / of the minimum ionization timescale from 
the rate equations. The minimum value found for At c is 
stored (as the new global time-step At). At this point we 
check which cells can be considered active during the next 
simulation step comparing the individual At c with At. If 
At c < ftsAt, where fts ^ 1 is a user-defined parameter, the 
cell is considered active and we assign At c = At. If a previ- 
ously active cell does not meet this criterion we de-activate 
the cell until a simulation time fts At (or a number nt s < fts 
of simulation steps) has elapsed. If nts = 1 all (leaf) cells 
are always considered active. 

Since the cells within the Ionization-front have the 
shortest At c , they are always selected as active. In prac- 
tice, the factor fts controls how broad is the region of 
active cells around the I-front. Cells outside of this re- 
gion are updated only when necessary, typically after sev- 
eral simulation time-steps, given their much longer ioniza- 
tion/recombination time-scales. Even for very large values of 
fts (e.g., fts = 10 4 ) the gain in computational speed is huge, 
since most of the volume has evolution time-scales several 
orders of magnitude larger than the I-fronts. 

3.3.3 Evolving the medium properties 

In this part of the computational algorithm, we solve the 
rate equation and we evolve the medium properties for the 
active cells with the newly calculated Ti and d . For the inte- 
gration of the rate equations, we use the Radau IIA method 
(Hairer & Wanner 1996), an implicit Runge-Kutta scheme 
of variable, adaptive order. This allows this part of the code 
to be computationally stable and reasonably fast for the re- 
quired accuracy. 

3.3.4 Adaptive refinement of the grid structure 

If activated, the Adaptive Mesh Refinement is performed 
with the recursive clustering algorithm of Berger & Rigout- 
sos (1991). Before applying this algorithm, we flag the cells 
that satisfy a given refinement criterion. 

The currently implemented refinement criterion is based 
on a combination of three variables for each atomic species: 
i) the cell neutral fraction, ii) the cell ionization rate (r*), 
and iii) the (frequency dependent) cell optical depth Ar c . 
The ultimate goal is to flag and refine the fast evolving cells 
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that are inside the I-front (i.e., the cells with the highest 
until their Ar c is below the desired value. In practice, the 
refinement criterion should produce a multi-mesh grid where 
the front is well resolved, i.e. where the Ar c at the species 
energy threshold is well below unity. Moreover, we would 
like also to have a pre-refined region just beyond the current 
position of the I-front, and a sharp cut corresponding to the 
regions where the I-front is just passed leaving the medium 
(highly) ionized. The size and shape of the refining region 
can be easily tuned with a criterion that includes the three 
variables discussed above, see Test 6 for an example. 

The recursive AMR algorithm starts from the base grid 
at the coarser level and it is composed by three main steps: i) 
flagging, ii) clustering, and iii) refming/unrefining. Once the 
cells to be refined at the coarser level (U — 0) are flagged 
as discussed above, they are clustered into new rectangu- 
lar patches. These are splitted until the clustering efficiency 
(i.e., the ratio between the number of flagged and total num- 
ber of cell in the new patch) is greater than a chosen mini- 
mum threshold, typically 80%. In order to avoid the creation 
of grids that are too small (and thus inefficient for the ray- 
tracing algorithm), it is also possible to chose a minimum 
size for the newly created patches. In most of the situations, 
we find that a minimum size of two (parent) cells gives the 
best balance between number of newly created cells and ef- 
ficiency of the algorithm. We then propagate the physical 
properties from the coarser grid (P) to the newly created 
patches at level P +1 (refinement). Currently, only straight 
injection refinement has been implemented. If a cell at level 
l* that was previously refined has not been flagged during the 
current time-step, the properties of its subcells at level 
are propagated back to the parent cell that becomes now a 
leaf cell (un-refinement). Temperature during un-refinement 
is weighted according to the subcell electron number den- 
sities, while the species neutral fractions (yi) are weighted 
according to the mass densities. This algorithm is recur- 
sively repeated, patch by patch, at finer levels, until there 
are newly flagged cells or the maximum number of refine- 
ment levels is reached. 

During clustering and refinement, the grid-tree and the 
cell-tree are updated, new memory is dynamically allocated 
for the newly created grids and, at the same time, the mem- 
ory slots associated with the un-refined grids are released. 
By the algorithm construction, parent and son grids have 
memory slots allocated in close parts of the memory map. 
This ensures an efficient memory usage and reduces compu- 
tational overhead during ray-tracing. 

In case where the AMR is performed on an existing 
multi-mesh grid (e.g., the output of an hydrodynamical 
AMR code), there is the very important option in RADAMESH 
to consider this initial grid as fixed: i.e., the AMR builds the 
new patches on the top of the existing multi-mesh grid with- 
out unrefining it below the original, initial level. This allows 
us to increase the resolution where needed for the radia- 
tive transfer (on the I-Front) without losing the resolution 
achieved with the hydro code. The capacity of RADAMESH of 
adding further levels of refinement on an existing AMR grid, 
without loosing previous information, is of great importance 
for several physical applications. For example for the study 
of the highly ionized (but denser than the average) regions 
in proximity of a QSO that where already touched by its 
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Figure 1. Test 1: image slices of the HI fraction at time t = 500 
Myr in the plane y = (top panel) and z = (bottom panel). 



I-front and that would appear as Lya forest in the quasar 
spectrum. 



4 VALIDATING TESTS 

In this section we present the validating tests of the code 
based on the radiative transfer code comparison project 
(Iliev et al. 2006; 106 thereafter). These tests have been de- 
signed in order to compare all the important aspects of sev- 
eral radiative-transfer codes present in the literature. They 
include the correct tracking of both slow and fast Ionization- 
fronts in homogeneous and inhomogeneous density fields, I- 
front trapping, spectrum hardening and the solution of the 
temperature state. The original tests in 106 (Test 1 to 4) are 
performed for a single, uniform grid, pure hydrogen medium 
and without recombination radiation (using the OTS ap- 
proximation). In order to compare our results with the other 
RT codes, we use the same single grid (pure- hydrogen) con- 
figuration as used in 106 to reproduce the results of Test 1 
to 4. In the second part of this section, we present a set of 
case studies aimed at substantiating the new characteristics 
of our code. In particular, we verify in Test 5 and Test 6 the 
multi-mesh, adaptive capability of RADAMESH. In Test 7 and 
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Figure 3. Test 1: time-evolution of the I-front position rj with 
respect to the Stromgen radius rg (bottom panel) and the ana- 
lytical solution (top panel). See text for details. 



Test 8 we show the effects of the diffuse radiation transfer 
for hydrogen only medium (Test 7) and including helium 
(Test 8). 

4.1 Test 1: isothermal HII region expansion 

This test represents the classical problem of the expansion of 
a Hff region in an uniform (pure-hydrogen) medium around 
a single ionizing source. We assume that a steady, monochro- 
matic ijnv — 13.6 eV) source emitting N 1 ionizing pho- 
tons per unit time turns on in an initially-neutral, uniform- 
density, static medium with hydrogen number density tih- 
Likewise in 106, we use for this test the OTS approximation. 
The temperature is fixed at T = 10 4 K. Under these con- 
ditions, and assuming that the front is sharp (i.e. that it is 
infinitely-thin, with the gas inside fully-ionized and the gas 
outside fully- neutral), there is a well-known analytical solu- 
tion for the evolution of the I-front radius, rr, and velocity, 
vi, given by 

1/3 



n 



VI = 



rs [1- 
rs 



- exp(-t/t rec ) 
exp (— t/t r , 



:) 



3trec [1 — exp( — t/t t 



12/3 



(16) 

(17) 
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Figure 5. Test 2: spherically averaged ionized (x) and neutral (1- 
x) fraction profiles at times t = 10 Myr (top panel) and t = 100 
Myr (bottom panel). The results from RADAMESH are presented 
as solid black lines. The other lines represent the results from a 
sample of four different codes taken from the RT code comparison 
project (106), in particular: C 2 -RAY (red, short-dashed line), CRASH 
(cyan, dotted line), FTTE (blue, long-dashed line), and RSPH (green, 
dot-dashed line). 



where 



rs = 



3N-, 



47ra B (T)n| 



1/3 



(18) 



is the Stromgren radius, i.e. the radius at which recombina- 
tions balance the ionizations and the HII region expansion 
stops. Here as(T) is the Case B recombination coefficient 
and 



[a fl (T)n H ] 



(19) 



is the recombination time. The HII region initially expands 
quickly and then slows down as the evolution time ap- 
proaches the recombination time, t ~ t lac . At a few recombi- 
nation times, the I-front stops and in absence of gas motions 
remains static thereafter. 

The numerical parameters for this test are the fol- 
lowings: computational box dimension L — 6.6 kpc (the 
source is at one corner of the box), gas number density 



Figure 6. Test 2: image slices of the gas temperature at coordi- 
nate z = and times t = 10 Myr (top panel) and t = 100 Myr 
(bottom panel). 



nji = 10 -3 cm -3 , initial ionization fraction (given by col- 
lisional equilibrium) x — 1.2 x 10~ 3 , and ionization rate 
N-y — 5 x 10 48 photons s . For these parameters the recom- 
bination time is t TCC = 3.86 x 10 15 s = 122.4 Myr. Assum- 
ing a recombination rate olb{T) — 2.59 x 10 _13 cm 3 s _1 at 
T = 10 4 K, then rs = 5.4 kpc. Note that the value of r s 
is actually independent from the use of the OTS or Case B 
approximation (see Test 7). 

In Figure [T] we show the images of the HI fraction in 
the plane y=0 and z=0 at time t = 500 Myr, when the 
equilibrium Stromgen sphere is reached. The HII region is 
nicely spherical in both the planes, demonstrating that our 
algorithm produces an uniform coverage of the solid angle 
around the source. In Figure [2j we plot the spherically av- 
eraged radial profiles of the ionized (a;) and neutral fraction 
(1-x) at times t — 30 and 500 Myr. The thickness of the 
transition between the HII and HI regions is in agreement 
with both analytical and numerical expectations from most 
of the other codes in 106. In particular, thanks to our new 
adaptive algorithm (that ensures the convergence of the ra- 
diation field in any cell), the I-front thickness does not suf- 
fer from diffusive effects shown by other "classical" Monte 
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Figure 7. Test 2: spherically averaged temperature profiles at 
times t = 10 Myr (top panel) and t = 100 Myr (bottom panel). 
The results from RADAMESH are presented as solid black lines. The 
other lines represent the results from a sample of four different 
codes taken from the RT code comparison project (106), see cap- 
tion in Figure [5] 

Carlo methods (cfr., e.g., the results of CRASH in 106). Fi- 
nally, in Figure |3J we show the time evolution of the I- front 
position (defined as the point of 50% ionization). The code 
tracks the I-front correctly, with the position never varying 
by more than few percent from the analytical solution. 

4.2 Test 2: HII region expansion with 
temperature evolution 

In this test we use the same parameters of Test 1, but now 
the ionizing source has a 10 5 K black-body spectrum and 
we allow the gas temperature to evolve. Initially, the gas is 
fully neutral with a temperature T=100 K. There are no 
analytical solutions for this test, therefore we compare our 
results to what obtained by the other codes in 106. 

In Figure[4] we show the images of the neutral hydrogen 
fraction (on the z — plane) at times t = 10 and t = 100 
Myr. The spherically averaged HI profiles, for the same time- 
snapshots, are presented in Figure [5] (black, solid lines) , to- 
gether with a sample of the results obtained by the other 
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Figure 8. Test 3: image slices of the HI fraction at coordinate 
z=0.5 (box units), at times t = 1 Myr (top panel) and t = 15 
Myr (bottom panel). 



codes in 106 (see caption in the Figure). The overall size of 
the HII region and the internal structure agree very well. 
The temperature images and the spherically averaged pro- 
files at times t = 10 and t = 100 Myr are presented in 
Figures [6] and [7] Also in this case, the resulting temperature 
structure agree well with most of the other codes in 106. In 
particular, we are able to obtain a large pre-heated region 
without any sign of spatial anisotropy, since our algorithm 
does not suffer from the under-sampling the radiation field 
(like in a "classical" Monte Carlo), even a large distances 
from the sources. 



4.3 Test 3: I-front trapping 

This test verifies that the propagation of an I-front within 
a dense clump is slowed down at the point of being stopped 
[trapped) and the production of a correct shadowing ef- 
fect (in both the ionization and the temperature state) be- 
hind the clump. We use the same set-up as described by 
106: a spherical (hydrogen-only) uniform cloud of radius 
r dump = 0.8 kpc, is located at the position (5, 3.3, 3.3) kpc 
within a box of length L — 6.6 kpc. The hydrogen num- 
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Figure 9. Test 3: image slices of the gas temperature at coordi- 
nate z=0.5 (box units), at times t = 1 Myr (top panel) and t = 15 
Myr (bottom panel). 



ber density and the initial temperature outside of the clump 



are, respectively, n out — 2 x 10 cm and T" u 



8000K, 



while inside the clump we have rii„ — 200n out = 0.04 cm 
and Tlfump — 40 K. The radiation has a black-body spec- 
trum with T e ff — 10 K and a constant ionizing photon flux 
F — 10 6 s -1 cm -2 , incident to the y = box side. 

Given these parameters, we expect that the I-front 
should be trapped slightly beyond the clump centre, as dis- 
cussed in 106. In Figure [8] we show slices of the gas neutral 
fraction at the box mid-plane z = 3.3 kpc (passing through 
the centre of the clump) at times t = 1 Myr (top panel) 
and £ = 15 Myr (bottom panel). The corresponding tem- 
perature slices are presented in Figure [9] At t = 1 Myr, the 
I-front is not yet trapped and it is still moving supersoni- 
cally from the edge of the clump. The images show that the 
shadow is sharp and produced correctly behind the clump. 
As expected, the I-front is trapped slightly beyond the clump 
centre at t = 15 Myr. The position of the I-front and the 
HI profiles inside the clump, presented in Figure |10| agree 
well with the other codes in 106, altough, especially at later 
times, the results are rather code dependent. The same is 



true for the temperature profiles, as shown in Figure (111 



Figure 10. Test 3: line cuts of the ionized and neutral fraction 
along the axis of symmetry through the centre of the clump at 
times t = 1 Myr (top panel) and t = 15 Myr (bottom panel). The 
results from RADAMESH are presented as solid black lines. The other 
lines represent the results from a sample of four different codes 
taken from the RT code comparison project (106), see caption in 
Figure [5] 



4.4 Test 4: Multiple sources in a "cosmological" 
density field 

In this test, we follow the propagation of the Ionization- 
fronts from multiple sources in a static cosmological density 
field. The initial conditions for the density and source spatial 
distribution/luminosities are provided by 106. In particular, 
we use a time-slice at z = 9 from a hydro-simulation with a 
box size of 0.5 h' 1 comoving Mpc and 128 3 cells. The initial 
temperature is fixed to T=100 K everywhere. The ionizing 
sources correspond to the 16 most massive halos in the box, 
with a luminosity proportional to the halo mass and a black- 
body spectrum with T = 10 K (see 106 for more details). 

In Figures [T2| and \l3\ we present slices of neutral hydro- 
gen fraction cut through the simulation box at coordinate 
z = 0.5, in box units, at time t = 0.05 Myr and t = 0.2 
Myr, respectively. Comparing our results with three of the 
four codes that performed this test in 106 (C 2 -Ray, CRASH 
and FTTE), we find a general agreement, although all the 
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Figure 11. Test 3: line cuts of the temperature along the axis 
of symmetry through the centre of the clump at times t = 1 
Myr (top panel) and t = 15 Myr (bottom panel). The results 
from RADAMESH are presented as solid black lines. The other lines 
represent the results from a sample of four different codes taken 
from the RT code comparison project (106), see caption in Figure 

m 



codes produce somewhat different morphologies. This gen- 
eral agreement is confirmed also examining the temperature 
slices in Figures [T4| and |T5| Here the differences between the 
codes presented in 106 is higher and due to different assump- 
tion about spectral hardening and partially to the different 
algorithm used. The spectral hardening effect, namely the 
increased temperature in the regions that are still mostly 
neutral, is traced in detail by RADAMESH. 




» 
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Figure 12. Test 4: image slices of the HI fraction at coordinate 
z=0.5 (box units) and time t = 0.05 Myr obtained by RADAMESH 
(top-left panel) and by three of the four codes that performed the 
same tests in 106, in particular: CRASH (bottom-left panel), C 2 -RAY 
(top-right panel), and FTTE (bottom-right panel). 




r 



Figure 13. Test 4: image slices of the HI fraction at coordinate 
z=0.5 (box units) and time t = 0.2 Myr obtained by RADAMESH 
(top-left panel) and by three of the four codes that performed the 
same tests in 106 (see caption in Figure [ill. 



4.5 Test 5: static multi-mesh 

With this test, not present in the original set of 106, we show 
and verify the multi-grid capability of RADAMESH in the case 
of a static multi-mesh. In particular we use the same ini- 
tial condition of Test 1 changing the original grid with three 
nested, concentric grids with 64 3 cells, corresponding to two 
levels of refinement (with a factor two and four increased 
resolution with respect to the base mesh). With this config- 
uration, the centre of the box has a resolution equivalent to 



a 256 3 cells grid. We limit the number of levels and the res- 
olution of the base grid in this test for illustrative purposes. 
The source position is (x — 0, y = 0, z — 0.5) in box units. 



In Figure (16 1, we show the image of HI fraction and the 



computational meshes corresponding to the slice at coordi- 
nate z = 0.5 and time t = 200 Myr. As we can see from the 
image, the front is tracked very well despite of the different 
grids (and resolution) , with no spurious effect introduced by 
the multi-grid structure. 
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Figure 14. Test 4: image slices of the gas temperature at co- 
ordinate z=0.5 (box units) and time t = 0.05 Myr obtained by 
RADAMESH (top-left panel) and by three of the four codes that per- 
formed the same tests in 106 (see caption in Figure |12[) , 



RADAMESH 




Figure 15. Test 4: image slices of the gas temperature at co- 
ordinate z=0.5 (box units) and time f = 0.2 Myr obtained by 
RADAMESH (top-left panel) and by three of the four codes that per- 
formed the same tests in 106 (see caption in Figure [l2"]l . 



4.6 Test 6: adaptively refining multi-mesh 

In the original set of tests in 106, the initial conditions (e.g., 
box size, hydrogen density, source luminosity) have been 
chosen in such a way to properly resolve the I-front with 
a single, uniform grid of 128 3 cells. For example, in the final 
test of 106, Test 4 (multiple sources in a "cosmological" den- 
sity field), the box size is fixed to 0.5 h~ comoving Mpc, cor- 
responding to ~ 70 physical kpc at the simulation redshift 
(z — 9) for h — 0.7. However, in a more realistic cosmologi- 
cal situation, e.g. for the study of the reionization process or 
the expansion of the I-front around a bright, high-redshift 
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Figure 16. Test 5: (Ionization-front expansion in a static multi- 
mesh) HI fraction image slice at time t = 200 Myr and coordinate 
z = 0.5 (box units). The multi-mesh structure in the plane z = 0.5 
is overlaid. 



QSO, the box size must be at least 2 orders of magnitude 
larger than Test 4 (see, e.g., Meiksin 2009). 

This test, not present in the original set of 106, demon- 
strate the ability of RADAMESH to resolve the I-front of a 
bright quasar in a large cosmological box, with the help of 
an adaptively evolving multi-mesh. In particular we show 
that properly resolving the I-front is essential to accurately 
predict the temperature state of the gas, both in the I-front 
and in the QSO HII region. For this purpose we run a series 
of simulation increasing the maximum level of refinement 
until convergence is reached for the temperature state of 
the gas. In the same spirit of Test 2, we use an uniform 
(hydrogen-only) medium with evolving temperature and we 
locate the source at one corner of the box. The box size is 
100 comoving Mpc, the medium density is equal to the mean 
hydrogen density of the universe and it evolves with redshift. 
The initial redshift is z = 7.5, the total evolution time is 10 8 
yr (corresponding to a final redshift z ~ 6.7). The source has 
a power-law spectrum with a (frequency) spectral slope of 
a v — 1.7, sampled in 60 logarithmically spaced bins from 1 
to 60 Rydberg. The total ionization rate is N = 10 57 ph s . 
These parameters are comparable to the respective quanti- 
ties associated with the observed high-redshift quasars in 
SDSS (e.g., Fan et al. 2006). 

For the assumed spectrum (which has a mean ionizing 
photon energy of ~ 2.4 Ryd) the photon mean-free path is 
A™ p ~ 6 physical kpc at the initial redshift for a fully neu- 
tral patch of gas composed of hydrogen only. To resolve the 
I-front properly we would ideally need that A™ p is sampled 
by at least two cells, i.e. a spatial resolution of ~ 3 physical 
kpc. With an uniform grid this would correspond to a 4096 3 
cells mesh, effectively out of the reach for current compu- 
tational facilities. Instead, we use a 128 3 base grid and 5 
additional levels of adaptive grid refinement that follow the 
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Figure 17. Test 6: (Ionization-front expansion in a evolving multi-mesh) HI fraction image slice at time t = 100 Myr and coordinate 
z = (box units). The inset shows the evolved multi-mesh (composed of 5 levels of refinement) that closely follows the Ionization-front. 



I-front expansion, effectively achieving the required spatial 
resolution. 

In Figure |17| we show an image slice of the hydrogen 
neutral fraction at t = 10 8 yr (in the quasar rest-frame). 
The I-front scale is so small compared to the box size that it 
visually appears perfectly thin in the large image. Zooming 
in by a large factor in a small region containing the I-front 
(see inset in the same figure), and now the I-front thickness 
appears. In box units, the I-front size (as measured from 
the points where the neutral fractions are 0.1 and 0.9) is 
roughly ~ 0.0025, i.e. ~ 32 physical kpc, corresponding to ~ 
5A'" 



mfp 



Overlaid, we show the adaptively refined multi-mesh 
structure that closely follows the I-front, with the highest 
level of refinement (I = 5) that encompasses the I-front itself. 
At this level of refinement, the I-front is fully resolved by at 
least 10 simulation cells. 

In Figure [18] we show how resolving the I-front changes 
the predicted temperature profile inside the QSO HII bubble 
and in the I-front region. From bottom to top, the lines rep- 
resent the temperature profiles obtained varying the maxi- 
mum level of refinement from lmax = (i.e., single uniform 
grid) to lmax = 5. As expected, temperature convergence 
is reached for l max = 4, i.e. at spatial resolution compara- 



ble or smaller than \™ p , indicating that we are effectively 
resolving the I-front. In fact, fully resolving the I-front im- 
proves the sampling of the spectral hardening (within the 
I-front itself) and a harder spectrum is able to produce a 
lar ge in crease in the gas temperature, as observed in Fig- 
ure IS ! The prediction of the correct temperature state in 
the I-front and in the QSO bubble is fundamental, e.g., for 
the study of the reionization epoch with the near-zone Lya 
forest (see e.g., Bolton et al. 2010) or with the I-front Lya 
emission (Cantalupo et al. 2008). 



3 Note that the current temperature of a cell within the HII 
region at t = 100 Myr is mainly determined at an earlier epoch, 
i.e., when the cell was located within the I-front. Therefore, also 
the temperature increase within the HII region can be ascribed 
to the (I-front) spectral hardening effect. Numerical effects on the 
gas temperature due to the different resolution are negligible, as 
we show in Appendix B. 
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where the first equality is only valid for a monochro- 
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Figure 18. Test 6: (Ionization-front expansion in a evolving 
multi-mesh) Temperature profiles around a bright source (a QSO) 
in a large, cosmological box (100 comoving Mpc) and I-front adap- 
tive mesh refinement. The QSO was turned on at z = 7.5 and 
keep it at a constant ionizing rate of N = 10 57 ph s — 1 for t = 10 8 
yr (the final redshift is thus z ~ 6.7). From bottom to top, the 
lines correspond to increasing maximum level of refinement (from 
Imax = to l ma x = 5). The base grid is composed by 128 3 cells. 
This plot clearly shows that the high-resolution achieved through 
AMR is essential to correctly resolve the I-front and, thus, to 
properly recover the gas temperature profile around the QSO in 
a cosmological context. 



4.7 Test 7: Diffuse ionizing emission from HII 
recombinations 

We repeat Test 1 dropping the OTS approximation used 
in 106 and including now the full radiative transfer of dif- 
fuse photons produced by hydrogen recombinations. Since 
all the cell in the Stromgren sphere may be now source of 
ionizing photons, we must change the original test config- 
uration, placing the (discrete) source at the centre of the 
box and increasing the box size by a factor of two. Apart 
from these modifications, we use the same parameter set as 
in Test 1. 

For a monochromatic spectrum, homogeneous medium 
and for a fixed temperature, we know that the radius of 
the Stromgren sphere (R a ) obtained with the OTS (Case 
B) approximation must be equal to the one obtained with 
the full radiative transfer of the diffuse radiation. This is a 
fundamental constraint deriving from photon conservation. 
Moreover, once the ionization equilibrium has been reached, 
there is an analytical approximation for the ratio between 
the diffuse and directional radiation field needed to compen- 
sate, at a given radius r, the recombinations, as found by 
Ritzerveld (2005): 



7 (r) _ r diff (r) 



jdh-( r ) 



T dir (r) 



1 - 



a;/o 



(20) 



matic spectrum. Note that eq. ( 20 1 is exact for an "outward- 
only" system, i.e. when the contribution to 7 dlff (r) is given 
by the recombinations at r' < r (indeed, 7 dlff (0) = in 
Ritzerveld's solution). Therefore, we expect that this ap- 
proximation should slightly underestimate 7 dlff for r <C R s 
(note, however that within this region the total radiation 
field will be dominated by 7 dlr ). 

In Figure |19[ we present the result obtained by 
RADAMESH including the full radiative transfer of diffuse pho- 
tons (black solid line). In very good agreement with the ana- 
lytical expectations (red dashed line), we found that the dif- 
fuse radiation field strength equals the directional field from 
the central source at r = 0.877? s , confirming that diffuse 
radiation is dominant in the outer parts of the Stromgren 
sphere. The full RT results produce, as expected, a slightly 
higher value of diffuse radiation in the central parts of the 
HII regions with respect to the "outward- only" solution. Be- 
cause of photon conservation, this is correctly balanced by a 
lower diffuse field at the edge of the Stromgren sphere when 
compared to the analytical approximation. 

In Figure |20| we show how the hydrogen neutral frac- 
tion profile (1 — x) at "equilibrium" (t = 500 Myr) is mod- 
ified by the full radiative transfer of diffuse photons (red 
solid line) with respect to the OTS (Case B) approximation 
(blue, long-dashed line) and the other extreme represented 
by the Case A approximation (black, short-dashed line). As 
expected, in the interior of the HII regions, where the gas 
is highly ionized, the mean-free-path of the diffusion pho- 
tons is very large and thus the profile obtained by the full 
RT closely follows the one obtained assuming Case A. How- 
ever, as the local opacity increases moving outward, diffuse 
photons start to be absorbed and the HI fraction tends to 
become closer to the OTS case. In the outer parts of the HII 
region (r > 0.877? s ), diffusion radiation becomes the domi- 
nant component and the gas is more ionized than in the OTS 
Case B approximation. However, in accordance with photon 
conservation, the final size of the Stromgren sphere is very 
close to the one obtained with Case B: the recombination 
radiation escaped from the inner region has been absorbed 
(creating an "excess") closer to the Stromgren radius, but 
the final photon balance is the same. 

As a remark, it is interesting to note that the OTS 
approximation is actually never valid inside the Stromgren 
sphere for this very simple, standard test (apart, obviously, 
at the Stromgren radius). Actually, Case A is a much better 
approximation if one is interested in the inner part of the HII 
region, or, analogously, in the early phases of the Ionization- 
front expansion. 



4.8 Test 8: Diffuse ionizing emission from HII, 
Hell and Helll recombinations 

The previous test has been limited to the very simple case of 
a monochromatic spectrum (for both diffuse and direct ra- 
diation) and a hydrogen only, fixed temperature medium. 
In this way, we were able to compare our results to the 
analytical prediction from photon conservation arguments. 
However, the real effect induced by the full RT of recombi- 
nation radiation is much more complex since its characteris- 
tic spectral energy distribution (SED) is in general different 
from the SED of the discrete sources. Moreover, the presence 
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Figure 19. Test 7: (diffuse emission from Hff) ratio between the 
ionization rate from diffuse HH recombinations and direct emis- 
sion (monochromatic spectrum) as a function of distance from the 
central source. For comparison, the analytical "outward-only" so- 
lution found by Ritzerveld (2005) is shown as a red dashed line. 
Note that the analytical solution should became exact in proxim- 
ity of R s (see text for details). In this test, the medium temper- 
ature has been fixed to T = 10 4 K 




Figure 20. Test 7: (diffuse emission from HII) HI fraction profile 
at t = 500Myr for Case A (short dashed black line), "on-the-spot" 
Case B (long dashed blue line) approximations and including the 
full radiative transfer of diffuse emission from HII recombinations 
(solid red line). In this test, the medium temperature has been 
fixed to T = 10 4 K. 



of helium may have a profound impact also on the neutral 
hydrogen distribution, given that both the bound-free and 
the Balmer continuum of Hel and Hell (together with two- 
photon and line emissions) are able to ionize HI. 

In this test, we verify these effects using the same config- 
uration of Test 2, including now the full RT of recombination 
radiation from hydrogen and helium. The only parameters, 
apart the inclusion of the diffuse component RT, that have 
been modified with respect to Test 2 are the following: i) the 
discrete source has been placed in the centre of the box and 




til 

0.2 0.4 0.6 

r (box units) 

Figure 21. Test 8: (diffuse emission from HII, Hell and Helll) 
HI fraction profile at i = 500 Myr for Case A (short dashed 
black line), Case B (long dashed blue line) approximations and 
including the full radiative transfer of ionizing diffuse emission 
from HII, Hell and Helll recombinations (solid red line). 



the box size has been increased by a factor 2, ii) the medium 
include hydrogen and helium, iii) we evolve the simulation 
to t — 500 Myr. Note that we do not expect in this case a 
full equilibrium to be reached, given the complex interplay 
between hydrogen and helium recombination emission and 
their ionization state. 

In Figure |2~l) we present the HI profile obtained by the 
full RT of diffuse radiation (red solid line) in comparison 
with the Case A approximation (black, short-dashed line) 
and the OTS Case B (blue, long-dashed line) at t = 500 
Myr. In Figure |22| we show the analogous profiles for the 
Hel, Hell and Helll fractions (see labels in the figure). As 
evident from Figure |21[ there is now an excess of HI pho- 
toionizations with respect to Case B, that leads to a larger 
Stromgren sphere size. This excess is due to diffuse pho- 
tons produced by helium recombinations, in particular the 
Balmer continuum (plus Lya and two-photon continuum) 
produced in the Helll region, that alters the total photon 
budget available for the photoionization of HI. In effect, 
a fraction of the high-energy, Hell-ionizing photons from 
the central source are "converted" via Helll recombinations 
to lower energy, Hi-ionizing photons (and also Hel-ionizing 
photons; indeed, a similar but less conspicuous effect is also 



visible for the Hel fractions in Figure 22 1. Therefore, as- 



suming Case B for hydrogen and helium actually results in 
the loss of this important component of the photon bud- 
get. This enlightens the importance of including Hell and 
Helll recombination emission also for studies that are only 
interested in the properties of the hydrogen component. 



5 DISCUSSION AND CONCLUSIONS 

We have presented a new, three-dimensional radiative trans- 
fer code, called RADAMESH, based on an adaptive Monte Carlo 
ray-tracing scheme. The algorithm has been specifically de- 
veloped to efficiently resolve the small scales associated with 
the Ionization-fronts within large, cosmological simulations, 



© RAS, MNRAS 000,[TJ{T9] 



RADAMESH 17 



0.1 



0.01 



r\HeI 


1 




, | , i i | 


" NX 








_ K. 

// 
// 


\\ 

X 




/ %\ 

S7 \\ - 


// 

z /, Hel 

: / 

r| 




\\ 
/ / 


1 
| 

! Hel/ 


/ / 
// 
// 

/ 


/ 




\ N X 


- // 
// 

- ,//, 


/ 




\\ 

A 

\i \ - 



0.2 0.4 0.6 

r (box units) 



Figure 22. Test 8: (diffuse emission from HII, Hell and Helll) 
Hel, Hell and Helll fraction profiles at t = 500 Myr for Case 
A (short dashed black lines), "on-the-spot" Case B (long dashed 
blue lines) approximations and including the full radiative trans- 
fer of ionizing diffuse emission from HII, Hell and Helll recom- 
binations (solid red lines). 



with the help of an Adaptive Mesh Refinement (AMR) 
scheme combined with a new "cell-by-cell Monte Carlo" ap- 
proach: rays are casted separately from selected cells within 
the computational box to the sources with a method that 
ensures their correct solid angle distribution. 

This method has several advantages with respect to a 
classical Monte Carlo ray-tracing scheme]^] in particular for 
multi-mesh domains. Indeed, in a classical Monte Carlo, the 
number of rays to cast is determined by the smallest, most 
distant cells to the source, i.e. from the (few) cells at the 
highest refinement level, oversampling the rest of the box. In 
RADAMESH, the number of rays is proportional to the number 
of cells to be evolved during the current time-step, i.e. the 
fast evolving cells typically located within the I-front. This 
translates into a huge gain in computational speed and effi- 
ciency, given the small fraction of the computational volume 
occupied by the I-front. Moreover, we are now able to choose 
a local criterion that adaptively determines the number of 
rays needed in a particular cell, e.g. set by the convergence 
of the radiation field. This ensures to obtain the required 
accuracy in the prediction of the temperature or ionization 
state of the gas in a very efficient way: concentrating the 
computational efforts only where actually needed. 

In the same spirit of hydrodynamical AMR codes, 
RADAMESH is also able to increase the spatial resolution 
where required, adaptively refining the mesh in correspon- 
dence of the Ionization-fronts. The implemented algorithm 
is based on the original patch-based AMR method of Berger 
& Rigoutsos (1991). 

RADAMESH is able to trace the Ionization-fronts from 
multiple sources as well as the diffuse ionizing radiation 
produced by hydrogen and helium recombinations with a 
multi-frequency approach. The time-evolution of six differ- 

4 where a series of rays is uniformly casted around a source to 
obtain the correct solid angle distribution. 



ent species (HI, HII, Hel, Hell, Helll, e) and the gas tem- 
perature is followed with a time-dependent, non-equilibrium 
chemistry solver based on an implicit Runge-Kutta scheme 
of variable, adaptive order. 

We performed all the four tests present in the radia- 
tive transfer code comparison project of Iliev et al. 2006, 
plus four additional, new tests aimed to substantiate and 
show the new characteristic of RADAMESH . The first four 
tests include the correct tracking of the I-front in homoge- 
neous (Test 1) and in- homogeneous density fields with mul- 
tiple sources (Test 4), I-front trapping behind a dense clump 
(Test 3), and the effect of photo-heating on the gas tempera- 
ture state (Test 2). These tests, like in the original 106 work, 
have been performed on a single, static mesh, without emis- 
sion from atomic recombinations in order to better compare 
our results with the other codes. 

RADAMESH results are in very good agreement with the 
majority of the other codes present in 106. In particular, 
thanks to the new adaptive algorithm, we have shown (Test 
1) that the recovered I-front thickness and structure does 
not suffer from any diffusive effects and numerical broad- 
ening typical of other Monte Carlo codes (given their poor 
sampling at large distances from the source). The same is 
true when the gas temperature is examined (Test 2). 

The second set of tests (Test 5 to Test 8) shows the 
ability of RADAMESH to deal with both static (Test 5) and 
adaptively evolving (Test 6) multi-mesh structures, and with 
diffuse radiation produced by hydrogen (Test 7) and helium 
(Test 8) recombinations. In Test 5, we reproduce the same 
results of the classical Stromgren sphere situation (Test 1) 
with a nested hierarchy of meshes at different refinement 
levels, showing that no spurious effects are introduced by 
the multi-mesh structure. 

Most importantly, we have shown in Test 6 that 
RADAMESH is able to fully resolve an expanding I-front from 
a bright source (e.g., a quasar) within a large, cosmologi- 
cal box (100 comoving Mpc size), thanks to an adaptively 
evolving mesh with several levels of refinement. The ability 
of resolving the I-front on cosmological scales is fundamen- 
tal for a large range of applications. For instance, recovering 
the correct gas temperature within the I-front and the QSO 
bubble is important, e.g., for the study of the reionization 
epoch with the near-zone Lya forest (see e.g., Bolton et al. 
2010) or with the Lya emission generated within the I-front 
(Cantalupo et al. 2008). In Test 6, we have demonstrated 
that a proper treatment of the spectral hardening inside a 
resolved I-front may result in a substantial increase of the 
gas temperature (AT ~ 10 4 K) within the whole HII region 
surrounding the bright, central source. 

The effect of diffuse radiation generated by hydrogen 
recombinations on the classic Stromgren sphere case (with 
monochromatic radiation) has been presented in Test 7. 
Here, we have verified that our treatment of diffuse radiation 
produces the same Stromgren sphere size with respect to the 
Case B approximation, in agreement with photon conserva- 
tion. Moreover, in accordance with the analytical prediction 
by Ritzerveld (2005), we have verified that the diffuse field 
becomes the dominant component at a distance from the 
source corresponding to 87% of the Stromgren radius. The 
neutral fraction profile in the inner part of the HII region 
follows the same profile obtained with Case A approxima- 
tion, while in the outer part the gas is more ionized and the 



© RAS, MNRAS 000,[T|fl9l 



18 S. Cantalupo & C. Porciani 



I-front narrower than the result obtained with the widely 
used, Case B approximation. In Test 8, we have also consid- 
ered the effect of Hel and Hell recombination radiation on 
the hydrogen and helium ionization state. We have found 
that Hell diffuse radiation, especially the Hi-ionizing, Hell 
Balmer continuum (together with Hell two-photon contin- 
uum and Lya emission) may have an important effect also 
on the hydrogen ionization state, substantially increasing 
the Stromgren sphere size. 

At present, RADAMESH is not yet coupled with hydro- 
dynamics but is able to post-process the output of three dif- 
ferent hydrodynamical codes, both grid-based, e.g., RAMSES 
(Teyssier 2002), CHARM (Miniati & Colella 2007) and particle- 
based, e.g. GADGET (Springel 2005). These outputs are effi- 
ciently converted to RADAMESH patch-based structures with 
a fast algorithm based on the clustering method of Berger & 
Rigoustous (1991), also used into the AMR module. The de- 
fault format is very similar to the widely used CH0MB0 HDF5 
structure^] allowing RADAMESH outputs to be directly visu- 
alized with the most recent, high-performance visualization 
packages (e.g., VISIlQl. Currently, RADAMESH is efficiently 
parallelized with OpenMP. All the results presented in this 
paper have been obtained with a few hours of computational 
time on a 8-core Intel-Xeon workstation. 
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APPENDIX A: ANALYTICAL 
APPROXIMATION FOR THE SOLID ANGLE 
OF A CUBIC CELL 

Starting from the explicit form of the solid angle for each face 
of a cubic cell, we have found the following approximation 
for the total cell solid angle as seen by any source inside (or 
outside) the computational box: 



© RAS, MNRAS 000, [T]fl9] 



RADAMESH 19 



= E E MAX 

n=l,2 5 3 £=-1,1 

+ v 



E 



(A in - S) 



x 



arctan 



i=l 



X ip + X iq + X in 



(Al) 



where e pqn is the three-dimensional Levi-Civita symbol (note 
the implicit sum over the indices p and q) , 
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is the 8x3 matrix containing the cell vertexes coordinates 
(translated to the Cartesian system with the source at the 
origin), xf u and x s ° urce are, respectively, the cell centre 
and the source coordinates, I is the linear cell size, Js,i is 
the 8x1 unit matrix (i.e., the matrix of ones), and A^.j is 
the following 8x3 matrix: 



-1 
-1 
1 
1 

-1 
-1 
1 
1 



All lengths and coordinates are in box units. Note that 
the presence of the factors containing Aij and e pqn reduce 
the total number of terms in the overall sum to four per cell 
face. The function MAX removes from the sum the solid an- 
gle of a face invisible to the source, since this has a negative 
sign. 

In Figure [A"T| we show the comparison between the re- 
covered cell solid angle obtained by a full Monte Carlo sim- 
ulation with 128 3 cubic cells and 10 s rays (black circles), 
and the above analytical approximation (red solid line). The 
solid angle is shown as a function of the distance from the 
source, placed at (1/2,1/2,1/2), and for two different direc- 
tions: parallel to the x-axis (i.e., normal to two cell faces; 
lower line) and along the line of sight that connects the 
source to the box vertex (1,1,1) (upper line). Note that 
these two directions correspond to the minimum and maxi- 
mum possible solid angle of a cubic cell at a given distance 
from the source. The analytical approximation is in very 
good agreement with the Monte Carlo results. 



APPENDIX B: TEST 6 FOR A SOURCE WITH 
SOFT SPECTRUM 

In Test 6, we have shown that the spectral hardening in- 
side a resolved I-front may result in a substantial increase 
of the gas temperature (AT ~ 10 4 K) within the whole HII 
region surrounding a bright quasar. To better demonstrate 
that this result is not affected by numerical effects due to the 
different grid resolutions, we have repeated Test 6 assuming 
a source with similar ionizing rate but with a soft spectrum 
(a black-body with T — 2 x 10 4 K). In this case, spectral 
hardening is minimal and therefore we do not expect that 
increasing the resolution on the I-front should substantially 
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Figure Al. Solid angle test: comparison between the recovered 
cell solid angle by a full Monte Carlo simulation (black circles) 
and our analytical solution (red solid lines) in a 128 3 grid. The 
solid angle is shown as a function of the distance from the source, 
as seen by two different angular directions (see labels). 
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Figure Bl. Same as Figure 18 (Test 6) but for a source with a 
soft spectrum (a black-body with T = 2 X 10 4 K). In this case, 
spectral hardening is minimal and different maximum levels of 
refinement produce now very similar temperatures inside the HII 
region. This demonstrates that numerical effects are not substan- 
tially affecting the results of Test 6. 



change the gas temperature as in Test 6. As we show in Fig- 
ure 



Bl (line colors and style have the same meaning as in 
Figure 18), this is indeed the case: different maximum levels 
of refinement produce now very similar temperatures inside 
the HII region showing that numerical effects are not impor- 
tant. Note that the (larger) shift in the I-front position in 
Figure Bl (with respect to Figure 18) is due to the fact that, 
at the lowest resolution, cells are very optically thick to ra- 
diation with frequencies right above the ionization threshold 
(the vast majority, given the soft spectrum assumed here). 
This slightly alters photon-conservation, which holds to high 
accuracy when the front is resolved into less optically thick 
elements, as shown in Test 1. 
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